Projected Density of States

This notebook will demonstrate use of the PDoS module for computing the Projected Density of States. We will show how it is used in both the cubic and linear scaling mode. After that, we’ll do some extra activities related to fragments. Finally, we will cover some current limitations of the O(N) mode.

geom = """HETATM    1  O   HOH     1      -0.702  -0.056   0.010  1.00  0.00           O
HETATM    2  H   HOH     1      -1.022   0.847  -0.011  1.00  0.00           H
HETATM    3  H   HOH     1       0.258   0.042   0.005  1.00  0.00           H
HETATM    4  O   HOH     2       2.221   0.027   0.001  1.00  0.00           O
HETATM    5  H   HOH     2       2.597  -0.412   0.767  1.00  0.00           H
HETATM    6  H   HOH     2       2.593  -0.449  -0.745  1.00  0.00           H  """

from io import StringIO
from BigDFT.IO import read_pdb
with StringIO(geom) as ifile:
    sys = read_pdb(ifile)

from BigDFT.Calculators import SystemCalculator
code = SystemCalculator(skip=True, verbose=False)

Cubic Scaling Mode

First, we will show how to plot the PDoS in the cubic scaling mode.

from BigDFT.Inputfiles import Inputfile
inp = Inputfile()
inp.set_hgrid(0.4)
inp.set_xc("PBE")

We need to activate the PDoS feature. We will also add some extra empty orbitals for later comparison with the linear mode, which automatically has some empty states.

inp.calculate_pdos()
inp.add_empty_scf_orbitals(4)

log = code.run(sys=sys, input=inp, name="pdos-C")

Once the calculation has been completed, we access which projections are available.

from BigDFT.PDoS import cubic_projection_info
proj_list = cubic_projection_info(log)
print(proj_list)

{'O': {'s': [0, 6], 'px': [1, 7], 'py': [2, 8], 'pz': [3, 9]}, 'H': {'s': [4, 5, 10, 11]}}

Then using this information, we can process the pdos.dat file to get the weights.

from BigDFT.PDoS import cubic_projection_weights
weights = cubic_projection_weights(log, proj_list)

Then it is time to create a DoS object to help with the plotting.

from BigDFT.DoS import DoS
dos = DoS(energies=log.evals[0][0], units='AU', label="Total")

from numpy import array  # Needed since we're hopping over the k points
dos.append(energies=array(log.evals[0]), units='AU',
           norm=[[x['H']['s'] for x in weights]], label="Hs")
dos.append(energies=array(log.evals[0]), units='AU',
           norm=[[x['O']['s'] for x in weights]], label="Os")
dos.append(energies=array(log.evals[0]), units='AU',
           norm=[[x['O']['px'] + x['O']['py'] + x['O']['pz'] for x in weights]], label="Op")

Add the sum as a sanity check.

dos.append(energies=array(log.evals[0]), units='AU',
           norm=[[x['H']['s'] +
                  x['O']['s'] + x['O']['px'] + x['O']['py'] + x['O']['pz']
                  for x in weights]], label="Sum")

Now we are ready to make a picture.

from bokeh.io import output_notebook
output_notebook(hide_banner=True)
def label_fig(p):
    p.add_layout(p.legend[0], 'right')
    p.xaxis.axis_label = "Energy (eV)"
    p.legend.click_policy = "hide"
    p.toolbar_location = "right"
    p.aspect_ratio = 2
    return p

from bokeh.plotting import figure, show
from bokeh.palettes import Category10

p = figure()
for i, k in enumerate(dos.labels):
    color = Category10[10][i % 10]
    p.line(dos.get_curves()[k].x, dos.get_curves()[k].y, line_color=color, legend_label=k)

show(label_fig(p))

Linear Scaling Mode

Now we can work on the linear scaling mode. First, the calculation has to be done. The support_function_multiples option is needed to help understand how well our support functions approximate s, p, etc.

inp = Inputfile()
inp.set_hgrid(0.4)
inp.set_xc("PBE")
inp["import"] = "linear"
inp["lin_general"] = {"support_function_multipoles": True}

log = code.run(sys=sys, input=inp, name="pdos-L")

To post-process the linear scaling mode, we have to diagonalize the matrix to get the eigenvectors and eigenvalues. We also need the density matrix for projecting.

from BigDFT.PostProcessing import BigDFTool
tool = BigDFTool()
hmat = tool.get_matrix_h(log)
smat = tool.get_matrix_s(log)
kmat = tool.get_matrix_k(log)

from scipy.linalg import eigh
evals, evec = eigh(hmat.todense(), b=smat.todense())

Now we can apply a similar processing as to the cubic mode to get the weights. Importantly, you want to pass the density (times one half) and overlap matrix in if you are using the Mulliken representation. Otherwise, the eigenvectors should be recomputed after Lowdin orthogonalization is applied.

from BigDFT.PDoS import linear_projection_info, linear_projection_weights
proj_info = linear_projection_info(log)
weights = linear_projection_weights(evec, proj_info, s = smat, k = 0.5 * kmat)

Let’s plot and see the results.

dosL = DoS(energies=evals, units='AU', label="Total")
dosL.append(energies=array([evals]), units='AU',
            norm=[[x['H']['s'] for x in weights]], label="Hs")
dosL.append(energies=array([evals]), units='AU',
            norm=[[x['O']['s'] for x in weights]], label="Os")
dosL.append(energies=array([evals]), units='AU',
            norm=[[x['O']['px'] + x['O']['py'] + x['O']['pz'] for x in weights]], label="Op")
dosL.append(energies=array([evals]), units='AU',
            norm=[[x['H']['s'] +
                   x['O']['s'] + x['O']['px'] + x['O']['py'] + x['O']['pz']
                   for x in weights]], label="Sum")

p = figure()
for i, k in enumerate(dosL.labels):
    color = Category10[10][i % 10]
    p.line(dosL.get_curves()[k].x, dosL.get_curves()[k].y, line_color=color, legend_label=k)
show(label_fig(p))

Notice how the density matrix has been used to filter out the virtual orbitals. They could be recovered by passing the overlap matrix as the density matrix.

We will also plot the linear and cubic values side by side to get a better idea of the differences.

p = figure()
for i, k in enumerate(dosL.labels):
    color = Category10[10][i % 10]
    p.line(dosL.get_curves()[k].x, dosL.get_curves()[k].y, line_color=color)
    p.line(dos.get_curves()[k].x, -1*dos.get_curves()[k].y, line_color=color, legend_label=k)
p.yaxis.axis_label = "Cubic <-- DoS --> Linear"
show(label_fig(p))